Effective-interaction approach to the many-boson problem 
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We show that the convergence behavior of the many-body numerical diagonalization scheme for 
strongly interacting bosons in a trap can be significantly improved by the Lee-Suzuki method adapted 
from nuclear physics: One can construct an effective interaction that acts in a space much smaller 
than the original Hilbert space. In particular for short-ranged forces and strong correlations, the 
method offers a good estimate of the energy and the excitation spectrum, at a computational cost 
several orders of magnitude smaller than that required by the standard method. 
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The many-body problem of interacting bosons or 
fermions poses a continuous challenge for quantum 
physics. The complexity of the quantum states increases 
quickly with higher particle numbers or stronger correla- 
tions, making it essentially impossible to solve the many- 
body problem exactly. Thus, mean-field methods are 
often applied, being numerically much less demanding 
than any attempt to diagonalize the many-body Hamil- 
tonian. Examples are the Gross-Pitaevskii approach for 
trapped bosons 0, 0, S|j or the celebrated Kohn-Sham 
equations 0, Q for fermions, as often used in the local 
(spin) density approximation. However, these methods 
can treat correlations only in an approximate way, and 
are proven insufficient for strong interactions between the 
particles. Quantum Monte Carlo calculations provide an 
alternative approach in many cases (see for example, the 
review by Harju 0), but have other drawbacks, such as 
limited accessibility of the excitation spectrum, a prob- 
lem shared by the other methods. The development of 
alternative diagonalization methods such as, e.g., coupled 
cluster methods Q is an urgent issue. 

Facing the above problems with strong interactions, 
one quickly realizes that often, the only method of choice 
is the straightforward numerical diagonalization of the 
many-body Hamiltonian. In fact, this has been tried for 
a large variety of physical problems, ranging from nu- 
clear structure (see for example, the review by Caurier 
et al. 0) and quantum chemistry to artificially made 
quantum systems such as metallic clusters [§] and quan- 
tum dots [10|| . However, with increasing particle number 
or stronger interactions between the particles, the num- 
ber of basis states needed for an accurate description of 
the many-body quantum system increases beyond com- 
putational reach. Truncations of Hilbert space become 
necessary - but often for the reduced basis, the results 
are too inaccurate [111 ]. 

In nuclear physics, a significant step forward has been 
achieved by applying the Lee-Suzuki method [H, 
that prescribes unitary transformations on operators (e.g. 
the Hamiltonian) to obtain effective operators within the 
reduced basis space. This method has been successfully 
used in so called no-core shell model calculations, for ex- 
ample, where nuclear systems with ~12 fermions have 



been studied |l4l |. 

Since the experimental realization of Bose-Einstein 
condensation with cold, trapped atoms, much interest 
turned to the physics of harmonically confined many- 
boson systems. Although the trapping potentials today 
confine thousands of bosons, the few-body regime does 
not seem to be impossible to reach, having in mind also 
the recent advances with optical lattices fla |. Increasing 
technological expertise with Bose-Einstein condensates 
on 

chips" (If}, 

together with the newly emerging research 
area of "atomtronics" fl7l |. makes the need for further the- 
oretical developments for a description of cold-atom gases 
with strong correlations an urgent issue. 

In this Letter, we apply the Lee-Suzuki method - to the 
best of our knowledge, for the first time - to a system of 
harmonically trapped spinless bosons with short-ranged 
repulsive interactions. The fact that the Lee-Suzuki 
method works very well for the strong short-ranged inter- 
actions between the nucleons encourages its application 
to describe cold atom gases beyond mean-field. 

At this point we note that there exist alternative 
ways of doing the unitary transformation. Two exam- 
ples, so called renormalization group transformations, are 
" Vin w-k" (l8| and SRG (similarity renormalizati on g roup) 
[191 ] - SRG is frequently used in nuclear physics [20|. An- 
other example, related to SRG, is UCOM (unitary corre- 
lation operator method) [2l[. In future studies, it would 
be interesting to try alternative methods for the prob- 
lem at hand, but here we have focused on the Lee-Suzuki 
transformation. Another study which applies the Lee- 
Suzuki transformation on a non-nuclear system is [22| . 
in which the two-electron problem in a quantum dot is 
examined. Also, in [H, 0] the few-body problem of a 
trapped Fermi gas is investigated with techniques similar 
to what is used here. 

For an ultra-cold gas of neutral atoms, the interaction 
between the particles can often be modeled by a short- 
range potential. We choose to use a Gaussian distribu- 
tion function parameterized by a range a, and a strength 
coefficient g. We here let the system be (quasi-) two- 
dimensional. The Hamiltonian is 
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Interaction strength Range Scattering length 

9 <x a 

1 0.1 0.000283 

10 0.1 0.0871 

10 1 0.871 

Table I: Numerically calculated s-wave scattering lengths for 
some parameter choices of the interaction potential (all quan- 
tities in oscillator units). 



1=1 j^Ll 

where TV is the number of particles and m is the parti- 
cle mass. The chosen interaction is normalized so that it 
becomes a (5-function in the limit when a goes to zero, 
which also allows us to avoid the mathematical difficulties 
of point-interactions in connection to the diagonalization 
of the full many-body Hamiltonian [25[ . In a cold atomic 
gas, the strength coefficient g is related to the scattering 
length between the particles. The scattering length can 
in some systems be experimentally tunable by Feshbach 
resonances 0]. The typical length scale of the system is 
the oscillator length, y/h/{muj). The range of the inter- 
action in ultra-cold atomic gases is typically very short 
compared to the w avelength of the particles; here we pick 
a = O.lyft/ (moj) for our calculations (in the following 
the oscillator length unit is used). This parameter will 
be discussed in more detail later. 

In two dimensions, the (s-wave) scattering length, a, is 
usually defined to be a positive quantitity, see for exam- 
ple Ref. jH]. A two-body system will have at least one 
bound state when the scattering length is larger than the 
distance outside which the interaction potential is zero. 
However, a purely repulsive short-range potential pro- 
duces a scattering length that is smaller than the range 
of the potential. For such systems, an increased attrac- 
tion (or reduced repulsion) would decrease a towards zero 
and an even stronger attraction would produce a bound 
state and give an a decreasing from infinity. The scat- 
tering lengths corresponding to the repulsive interaction 
potential used here can be calculated numerically, see ta- 
ble m 

The short-range interaction between the particles in- 
duces short-range correlations, thus one would need a 
very large set of basis states in order to accurately de- 
scribe the wavefunction and to obtain a good estimate of 
the energy of the system. However, one may perform a 
unitary transformation of the Hamiltonian so that states 
within a given model space become decoupled from the 
ones in the complementary excluded space [13, ■ The 
finite model space and the complementary infinite space 
can be defined by the projection operators P and Q, re- 
spectively. The original Hamiltonian of the system is 
transformed to an effective Hamiltonian acting only in a 
subspace of the complete Hilbert space, while preserving 



a subset of the original eigenvalues. Thus, given the effec- 
tive Hamiltonian, the computational effort in finding the 
eigenvalues can be reduced. However, even though the 
original Hamiltonian contains only one- and two-body 
terms, the effective Hamiltonian will in general be an 
iV-body operator. So far, no approximation has been 
introduced, only an operator transformation given by: 

ft _ P + P$Q . QjP + P 

^ eft — i ^ original , j V / 

where £ is an operator acting as a mapping between 
the P- and Q-spaces, satisfying £ = Q£,P. Note that 
the transformation (HJ is unitary and yields an effective 
Hamiltonian that is energy-independent and hermitian. 
H e s does not couple states in the P-space with states in 
the Q-space. 

Unfortunately, to find the mapping operator £ one 
needs the exact solution of the iV-body problem, which 
is in itself the final goal. Our simplest, yet nontrivial, 
approximation is to develop a two-body effective Hamil- 
tonian. The approximation consists in finding £2 for the 
two-body problem and to compute an effective Hamilto- 
nian for the two-body system. By subtracting the one- 
body terms (kinetic + potential energy) one thus obtains 
an effective two-body interaction. This effective inter- 
action is then used to construct the TV-body effective 
Hamiltonian that will now contain only one- and two- 
body terms. Due to the approximation of the effective 
Hamiltonian, the obtained energies will not be bound 
by the variational theorem. This particular property of 
the current method is in contrast to the situation en- 
countered when using standard configuration interaction 
calculations, or quantum Monte-Carlo approaches. 

Although this approach can lead to good estimates of 
the energy eigenvalues, the obtained eigenvectors are in- 
correct since they do not contain components from the 
excluded Q-space. However, the eigenvectors are of- 
ten not interesting by themselves; only observables, ex- 
pressed as expectation values of physical operators, are 
relevant. Transformations similar to the one performed 
on the Hamiltonian, can be applied to other operators 
so that their expectation values can be obtained from 
the P-space eigenvectors. In the present study, however, 
we restrict the discussion to the many-body energy spec- 
trum. 

For bosons, the many-body basis states are perma- 
nents, with N particles distributed over the single- 
particle orbitals of the system. The truncation of the 
infinite basis is performed in the following way: Neglect- 
ing interactions, the state of lowest energy is the one 
with all N bosons in the lowest orbital as defined by the 
confinement potential, with energy E = HuN. We in- 
corporate in the P-space all many-body states with an 
energy 

E <hu(N + A/ max ), 
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so that Amax is a parameter determining the maximum 
allowed energy of particle-hole excitations from the state 
of lowest energy. The number of included states increases 
rapidly with A/" max . 

In addition, all possible combinations of two particles 
found in the A^-body model space (P) defines the re- 
stricted space of the two-body system (P2) that is used to 
compute the effective interaction. Having defined P2, the 
Hamiltonian for the complete two-body space is trans- 
formed. In practice, the complete two-body space must 
of course also be truncated. To this aim we take all two- 
body states that can be constructed using the first 20 
harmonic oscillator shells. 

The original Hamiltonian preserves angular momen- 
tum (L), implying that we can restrict the basis to states 
with a given value of L, thus limiting its size. We con- 
firm that the effective interaction also preserves angular 
momentum (within numerical accuracy). 

First we use the method to calculate energies of a sys- 
tem with only four particles, N = 4, with an interaction 
range a = 0.1. Here, we consider only states with zero 
angular momentum. For such a small system, the stan- 
dard configuration interaction method can be applied. 
Figure Q] shows results from both types of calculations. 
For g — 1, both methods give good energy estimates 
already at small A/ max . But for g = 10 the standard cal- 
culations require a large basis set to obtain a reasonably 
converged energy estimate, while the effective interaction 
provides an answer with significantly less computational 
effort. Typically, for an A/" max around 5, the basis size is 
of order 10, while for A/" max around 25 it is of order 10 4 . 
These numbers depend strongly on the system parame- 
ters N and L, though. 

Note that for J\f max = 0, the standard configuration 
interaction calculation reduces to a perturbative calcula- 
tion, where the energy grows linearly with the strength of 
the interaction, g. In addition, for sufficiently large A/" max 
the two methods will be equivalent, since the P-space is 
then actually the full space. 

The method works well both for the ground state en- 
ergy and the excitations. An interesting observation is 
that for g = 10, at A/" max = 6 the standard calculation 
shows a degeneracy in the first excited state. This effect 
is spurious, however, since the energies split up when a 
larger basis is used. When using the effective interaction, 
this false degeneracy does not occur. 

Let us now turn to systems with larger number of par- 
ticles. Figure [2] shows calculations for TV = 9 bosons. In 
the case of g = 1 the effective interaction method again 
very rapidly produces an accurate estimate of the en- 
ergy. For g = 10, where the correlations in the systems 
are much stronger than for g = 1, the energies obtained 
when using effective interactions appear to have reached 
some plateau, but still show a slow decrease with grow- 
ing A/" max - In comparison, the energies obtained from the 
standard calculations show a much slower convergence. 

The convergence behavior for an even larger system of 
N = 20 particles is shown in Fig. 03 In the case of g = 1 



N = 4 particles, L = Oh, <7 = 0.1 
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5 10 15 20 25 5 10 15 20 25 
Nmax (Hilbert space cutoff) 

Figure 1: Energies for a system of N = 4 bosons and total 
angular momentum L = Oh, for different cutoffs of the many- 
body Hilbert space (parameterized by A/" ma x). The range of 
the interaction is a = 0.1, and two different strengths (g) 
are shown. The blue dashed curves are the results from stan- 
dard configuration interaction calculations, while the red solid 
curves are energies obtained using the effective interaction ap- 
proach. While the standard calculations require a large basis 
set to give a good energy estimate, the effective interaction 
provides roughly the same answer with significantly less com- 
putational effort. (Please note that for small A/" ma x, there are 
very few states in the P-space, and consequently only a few 
eigenvalues can be obtained.) 

N = 9 particles, L = Oh, CT= 0.1 




10 15 20 5 10 
Nmax (Hilbert space cutoff) 



15 20 



Figure 2: Same as Fig. [T] but for N = 9 bosons. Note the 
rapid convergence of energies using the effective interaction 
approach. For example, we note that in this case a calculation 
with A/"max = 6 requires 12 basis states, while A/" ma x = 20 
would correspond to 81097 states. 



it is possible to reach fairly converged results with both 
methods. When g = 10 neither method is able to provide 
converged energies for the range of A/" max considered here. 
An intuitive interpretation would be that A/" max should 
be of the same order as N, allowing about one extra unit 
of energy per particle, so that every particle has some 
freedom to adjust in the system. Since the number of 
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N = 20 particles, L = Oh. (7= 0.1 
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Figure 3: Same as Fig Q] but for N — 20 bosons. For g — 10, 
the energy obtained using effective interaction cannot be said 
to be converged as a function of jVmax. 



basis states grows very rapidly with both N and A/" max , 
it seems that the method is most suitable in the strongly 
correlated regime; and makes it possible to study larger, 
but not much larger, systems than with the standard 
method. 

All results shown so far involved the L = OH (non- 
rotating) states of the different systems studied. How- 
ever, the method is not restricted to the non-rotating 
case. Figure S^a) shows the energies for a system of 
N = 9 bosons, at angular momentum L = 9ft.. Here, 
the energy is not a strictly decreasing function of Afmax, 
but aside from this the convergence is similar to that seen 
in figure [21 An energy obtained with the standard config- 
uration interaction method is always an upper bound to 
the true value, since it is a variational approach. As the 
size of the basis space is increased, a better or equally 
good estimate is found. Calculations using an approx- 
imated effective Hamiltonian, as in this study, are not 
variational so higher-order terms may contribute with ei- 
ther sign to the energy. 

The Lee-Suzuki approach was invented to handle the 
short-range correlations in nuclei [H, US], and the short 
range of the interaction is known to be essential for the 
performance of the method. For all results presented this 
far, we have set the range parameter a = 0.1. Figure|D[b) 
shows results with a larger range, a = 1, for a system 
with N = 9 particles. Apart from some deviations for 
small A^max, the effective interaction here does not give 
an improved convergence rate compared to the standard 
calculations. This result suggests that for a smaller than 
0.1 the method would perform even better than demon- 
strated in e.g. figure [21 However, the numerical solution 
of the two-body problem becomes more difficult with de- 
creasing a, and our present implementation prevents us 
from exploring smaller a. 

As mentioned, the correct effective Hamiltonian would 



in general be an A-body operator. It would be interesting 
to examine how inclusion of e.g. an effective three-body 
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Figure 4: (a) Energy levels (in units of huj) for a system of 
N = 9 bosons and angular momentum L = 9ft (in this case 
there are no possible basis states for Amax < 9). (b) Energies 
for a system with N = 9 and interaction range a = 1 (L = Oft 
states). See caption of figure [T] for explanations. 



interaction would affect the numerical convergence, al- 
though the two-body approximation is exact in the limit 
N m ax — > co. As seen in our results, the obtained energies 
typically converge rapidly as functions of N max , imply- 
ing that an effective two-body Hamiltonian is sufficient 
in many situations. 

To summarize, in order to calculate properties of cold, 
bosonic atomic gases with strong short-range interactions 
between the particles, we have employed a unitary trans- 
formation of the Hamiltonian. The transformation is 
used to recover correlation effects which would otherwise 
be lost within the heavily truncated basis space we con- 
sider. The method is in practice a modification of the 
standard configuration interaction approach. In many 
cases it produces a good estimate of the energy, with a 
computational effort which is several orders of magnitude 
smaller than that required by the standard method. The 
main advantage, compared to many other approaches, is 
the accessibility of the excitation spectrum with signifi- 
cantly reduced numerical effort. 
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